//g++ getMCtrueFromMOM.C  -g -O3 -Wall -Werror -I$ROOTSYS/include -I$ALICE_BUILD/include -I$ALICE_ROOT/TPC -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TRD/Cal  -L$ROOTSYS/lib -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic -lMinuit -lEG -lGeom -lVMC -lProof -lProofPlayer -lXMLParser -lXMLIO -lSpectrum -lTreePlayer -lMLP -lGui -L$ALICE_BUILD/lib/tgt_linuxx8664gcc -lCDB -lSTEER -lRAWDatarec -lESD -lSTEERBase -lANALYSIS -lRAWDatabase  -lANALYSISalice -lAOD -lTPCrec -lTPCbase -lTRDbase -lTRDrec -lSTAT  -lASImage  -o getMCtrueFromMOM

//getMCtrueFromMOM file ktail
//fit with guas (or with tail if ktail ==1)
//so, e.g.: 
// ./getMCtrueFromMOM histRec-anchor10h_MC_PbPbLHC11a10b_bis-091120122110hMChighpt_RAA_LARGE_SECO.root 1 > see1.log
// ./getMCtrueFromMOM histRec-anchor10h_MC_PbPbLHC11a10b_bis-091120122110hMChighpt_RAA_LARGE_SECO.root 0 > see0.log

/*
#include <stdio.h>

#include "TTreeStream.h"
#include "TMath.h"
#include "TH2D.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"

const Double_t EPSILON = 1e-12;
*/

#include "AnaDef.h"
#include "AnaAliDef.h"
#include "XGLUtils.h"
#include "dEdxGausN.h"

Double_t ResolutionGaus0(const Double_t xx, const Double_t amp, const Double_t mean, const Double_t resolution)
{
  //npar = 3
  return amp*TMath::Gaus(xx, mean, mean*resolution, kTRUE);
}

Double_t ResolutionGaus(const Double_t xx[], const Double_t par[])
{
  return ResolutionGaus0(xx[0], par[0], par[1], par[2]);
}

Double_t TailedGaus0(const Double_t xx, const Double_t amp, const Double_t mean, const Double_t resolution, const Double_t lambda)
{
  //
  //exp(-x/lambda)
  //
  
  const Double_t sigma = resolution*mean;

  Double_t sum=0, sW=0;

  //lambda=0, exp(-x/lambda)/lambda -> delta function
  if(lambda<EPSILON){
    return amp*TMath::Gaus( xx, mean, sigma, kTRUE);
  }

  //nlam < 5 step not enough to cover to tail to 1e-4
  //nlam < 10 not enough to 1e-5
  for(Double_t nlam=0; nlam<15; nlam+=0.5){
    const Double_t weight = TMath::Exp(-nlam);
    sum += TMath::Gaus( xx - nlam*lambda, mean, sigma, kTRUE)*weight;
    sW+=weight;
  }

  if (sW>EPSILON) return amp*sum/sW;
  return 0;
}

Double_t TailedGaus(const Double_t xx[], const Double_t par[])
{
  return TailedGaus0(xx[0], par[0], par[1], par[2], par[3]);
}

const Double_t mass[]={0.938272, 0.139570, 0.493677, 5.109989e-04 };
const char* tname[]={"proton","pion","kaon","electron"};
const char* vname[]={"fraction","mean","resolution","lambda"};
const char* mname[]={"integral","mean","rms","kurtosis"};
//test->
const Int_t ntype = 4;
//test<-

int main(int argc, char * argv[])
{
  //===================================================================================
  //                                      initialization 
  //===================================================================================

  //---- input arguments

  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc<3){
    printf("argc<3\n");return 1;
  }

  const TString fin(argv[1]);
  TFile::Open(fin);

  const Bool_t ktail(atoi(argv[2]));
  Bool_t kp2pt = kFALSE;
  if(argc==4)
    kp2pt=atoi(argv[3]);

  //---- controls

  dEdxGausN::HistPrint();

  //---- input hists

  TList *ll=new TList;
  
  TH2D *hrawdata=(TH2D*)gDirectory->Get("horiginal");
  if(!hrawdata){
    printf("no hh!\n");gDirectory->ls();exit(1);
  }
  /*
  //test-->
  hrawdata->Scale(0);
  TRandom3 ran3(1);
  //test<--
  */

  TH1D *hrawyield = (TH1D*)gDirectory->Get("horiginalyield");
  if(!hrawyield){
    printf("no hhrawyield!\n");gDirectory->ls();exit(1);
  }
  TH1D *hmeanx = (TH1D*)gDirectory->Get("meanx");
  if(!hmeanx){
    printf("no hhmeanx!\n");gDirectory->ls();exit(1);
  }
  ll->Add(hrawdata);
  ll->Add(hrawyield);
  ll->Add(hmeanx);

  TH2D * hpars[4]={0x0,0x0,0x0,0x0};
  for(Int_t itype=0; itype<ntype; itype++){
    //mcpar0;1        mcpar0 tpcsig : log10p (inequal x bin width)
    hpars[itype]=(TH2D*)gDirectory->Get(Form("mcpar%d",itype));
    if(!hpars[itype]){
      printf("no hpars %d\n",itype);gDirectory->ls();exit(1);
    }
    ll->Add(hpars[itype]);

    hpars[itype]->RebinY(2);
    printf("binwidth %f\n", hpars[itype]->GetYaxis()->GetBinWidth(1));
  }

  //---- output

  TTreeSRedirector *fStream = new TTreeSRedirector(Form("MCtrue_save_ktail%d.root",ktail));
  
  Int_t nfit = 0;

  //in case one of the 4 types is not updated due to small statistics like electrons, the value from last ix is used
  //if all types are not updated, the ix is just skipped    
  Double_t ys[4]={0,0,0,0};
  Double_t p0s[4]={0,0,0,0};
  Double_t p1s[4]={0,0,0,0};
  Double_t p2s[4]={0,0,0,0};
  Double_t p3s[4]={0,0,0,0};
    
  //===================================================================================
  //                                      Loops
  //===================================================================================

  //---- for each xx slice
  
  for(Int_t ix=1; ix<=hpars[0]->GetNbinsX(); ix++){
    //test-->
    //if(nfit>0) break;
    //test<--
    
    Double_t binx = hpars[0]->GetXaxis()->GetBinCenter(ix);
    Double_t meanx = hmeanx->GetBinContent(ix);
    if(meanx<dEdxGausN::Xmin())
      continue;

    if(dEdxGausN::FilterIX(ix))
      continue;

    Double_t binw = hmeanx->GetBinWidth(ix);
   
    nfit++;
    Double_t e0s[4]={0,0,0,0};
    Double_t e1s[4]={0,0,0,0};
    Double_t e2s[4]={0,0,0,0};
    Double_t e3s[4]={0,0,0,0};

    //totn only for 4 particles
    Double_t totn = 0;
    Double_t chi2[4]={0,0,0,0};
    Double_t ndf[4]={0,0,0,0};
      
    Double_t m0s[4]={0,0,0,0};
    Double_t m1s[4]={0,0,0,0};
    Double_t m2s[4]={0,0,0,0};
    Double_t m3s[4]={0,0,0,0};

    //---- for each type

    for(Int_t itype=0; itype<ntype; itype++){
      
      //---- get dE/dx distribution

      TH1D *h1d =hpars[itype]->ProjectionY(Form("pyx%dt%d",ix,itype),ix, ix);

      m0s[itype] = h1d->Integral(1,h1d->GetNbinsX());
      m1s[itype] = h1d->GetMean();
      m2s[itype] = h1d->GetRMS();
      m3s[itype] = h1d->GetKurtosis();

      if(m0s[itype]<1){
        continue;
      }
      ll->Add(h1d);

      //---- fit dE/dx 

      TF1 * ftmp=0x0;

      if(ktail)
        ftmp = new TF1(Form("f%d%s",ix,tname[itype]), TailedGaus,0,200,4);
      else
        ftmp = new TF1(Form("f%d%s",ix,tname[itype]), ResolutionGaus,0,200,3);

      ll->Add(ftmp);
      
      ftmp->SetParameters(h1d->Integral(0,1000), h1d->GetMean(), h1d->GetRMS()/ h1d->GetMean(), h1d->GetRMS()/10.);

      if(!kp2pt){
        h1d->Fit(ftmp,"L");
        printf("fitting done...\n\n");

        printf("fitresult %f %f --  %f %f %f %f\n", hrawyield->GetBinContent(ix),  h1d->GetXaxis()->GetBinWidth(1),  ftmp->GetParameter(0), ftmp->GetParameter(1), ftmp->GetParameter(2), ftmp->GetParameter(3));
      }

      //p0: counts
      p0s[itype]=ftmp->GetParameter(0)/h1d->GetXaxis()->GetBinWidth(1);
      p1s[itype]=ftmp->GetParameter(1);
      p2s[itype]=ftmp->GetParameter(2);
      p3s[itype]=ftmp->GetParameter(3);

      e0s[itype]=(kp2pt?sqrt(p0s[itype]): ftmp->GetParError(0))/h1d->GetXaxis()->GetBinWidth(1);
      e1s[itype]=ftmp->GetParError(1);
      e2s[itype]=ftmp->GetParError(2);
      e3s[itype]=ftmp->GetParError(3);

      totn += p0s[itype];

      /*
      //test->
      for(Int_t iran=0; iran<p0s[itype]; iran++){
        hrawdata->Fill(meanp, ran3.Gaus(p1s[itype], p1s[itype]*p2s[itype]));
      }
      //test<-
      */

      //---- get chi2 and ndf 

      for(Int_t ixtmp=1; ixtmp<=h1d->GetNbinsX(); ixtmp++){
        const Double_t tmpx = h1d->GetBinCenter(ixtmp);
        const Double_t yy = ftmp->Eval(tmpx);
        const Double_t cc = h1d->GetBinContent(ixtmp);
        if(cc<EPSILON)
          continue;

        const Double_t ee = h1d->GetBinError(ixtmp);
        chi2[itype] += pow(yy-cc,2)/ee/ee;
        ndf[itype]++;
      }

      ndf[itype]-= ftmp->GetNpar();

      //---- get hfit

      TH1D * hfit=new TH1D(Form("hfit%d%s",ix,tname[itype]), "", 500, h1d->GetXaxis()->GetXmin(), h1d->GetXaxis()->GetXmax());
      for(Int_t ixtmp=1; ixtmp<=hfit->GetNbinsX(); ixtmp++){
        const Double_t tmpx = hfit->GetBinCenter(ixtmp);
        const Double_t yy = ftmp->Eval(tmpx);
        hfit->SetBinContent(ixtmp, yy);
      }
      ll->Add(hfit);
    }
    if(totn<dEdxGausN::PbinThres()){
      continue;
    }

    //---- converting

    //p0: count -> fraction
    for(Int_t itype=0; itype<ntype; itype++){
      ys[itype]=p0s[itype];
      p0s[itype] = p0s[itype]/totn;
      e0s[itype] = e0s[itype]/totn;
    }

    //raw: internal fith, mn2d parameters
    Double_t raw0s[]={-999,-999,-999,-999};
    Double_t raw1s[]={-999,-999,-999,-999};
    Double_t raw2s[]={-999,-999,-999,-999};
 
    Int_t idx[4];
    TMath::Sort(ntype, p1s, idx, kFALSE);
    for(Int_t ith = 0; ith<ntype; ith++){
      raw0s[ith] = p0s[idx[ith]];
      raw1s[ith] = p1s[idx[ith]];
      raw2s[ith] = p2s[idx[ith]];

      if(ith>0){
        raw1s[ith] -= p1s[idx[ith-1]];
      }
    }

    //---- saving

    Double_t meanp = pow(10,meanx);
    Double_t binp = pow(10, binx);
    Double_t meanbg[4];
    Double_t binbg[4];

    for(Int_t itype=0; itype<ntype; itype++){

      meanbg[itype] = meanp/mass[itype];
      binbg[itype]  = binp/mass[itype];

      printf("check %d %d %e %e -- %e %e %e\n", ix, itype, p0s[itype], totn, p1s[itype], p2s[itype], p3s[itype]);
      
      (*fStream)<<"mn2d"<<
        Form("y%d=",itype)<<ys[itype]<<
        Form("rawp%d0=", itype)<<raw0s[itype]<<
        Form("rawp%d1=", itype)<<raw1s[itype]<<
        Form("rawp%d2=", itype)<<raw2s[itype]<<
        Form("p%d0=", itype)<<p0s[itype]<<
        Form("le%d0=",itype)<<e0s[itype]<<
        Form("he%d0=",itype)<<e0s[itype]<<
        Form("p%d1=", itype)<<p1s[itype]<<
        Form("le%d1=",itype)<<e1s[itype]<<
        Form("he%d1=",itype)<<e1s[itype]<<
        Form("p%d2=", itype)<<p2s[itype]<<
        Form("le%d2=",itype)<<e2s[itype]<<
        Form("he%d2=",itype)<<e2s[itype];
        
      (*fStream)<<"tree"<<
        Form("chi2%s=", tname[itype])<<chi2[itype]<<
        Form("ndf%s=", tname[itype])<<ndf[itype]<<
        Form("p%s%s=", tname[itype], vname[0])<<p0s[itype]<<
        Form("e%s%s=", tname[itype], vname[0])<<e0s[itype]<<
        Form("p%s%s=", tname[itype], vname[1])<<p1s[itype]<<
        Form("e%s%s=", tname[itype], vname[1])<<e1s[itype]<<
        Form("p%s%s=", tname[itype], vname[2])<<p2s[itype]<<
        Form("e%s%s=", tname[itype], vname[2])<<e2s[itype]<<
        Form("p%s%s=", tname[itype], vname[3])<<p3s[itype]<<
        Form("e%s%s=", tname[itype], vname[3])<<e3s[itype]<<

        Form("%s%s=", mname[0], tname[itype])<<m0s[itype]<<
        Form("%s%s=", mname[1], tname[itype])<<m1s[itype]<<
        Form("%s%s=", mname[2], tname[itype])<<m2s[itype]<<
        Form("%s%s=", mname[3], tname[itype])<<m3s[itype]<<

        Form("meanbg%s=", tname[itype])<<meanbg[itype]<<
        Form("binbg%s=", tname[itype])<<binbg[itype];  
    }
    (*fStream)<<"mn2d"<<
      "ix="<<ix<<
      "meanx="<<meanx<<
      "binx="<<binx<<
      "xx="<<binx<<
      "totn="<<totn<<
      "\n";

    (*fStream)<<"tree"<<
      "ix="<<ix<<
      "meanp="<<meanp<<
      "binp="<<binp<<
      "binw="<<binw<<
      "totn="<<totn<<
      "\n";
  }

  ll->ls();

  TFile *fout=fStream->GetFile();
  fout->cd();
  ll->Write();

  delete fStream;
}



